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ABSTRACT 

A theoretical method for determining equilibrium interface 
configurations in axisymmetric containers of arbitrary shape was 
investigated. A single differential equation was derived from the 
principle of minimum surface and potential energy using the calculus 
of variations. This equation, in conjunction with boundary conditions 
dependent on container shape and contact angle, can be numerically 
solved for the desired surface profile using the Runge-Kutta iteration 
technique. The method imposes no significant limitations on contact 
angle or Bond number and is easily programmed for computer solutions. 
Representative theoretical results are presented concerning the influ- 
ences of contact angle, Bond number, and container fill level on sur- 
face shapes. Also, theoretical results are compared with experimental 
data. 


NASA - GEORGE C. MARSHALL SPACE FLIGHT CENTER 



NASA-GEORGE C. MARSHALL SPACE FLIGHT CENTER 


TECHNICAL MEMORANDUM X- 53790 


LOW GRAVITY LIQUID- VAPOR INTERFACE 
SHAPES IN AXISYMMETRIC CONTAINERS 
AND A COMPUTER SOLUTION 

By 

Leon J. Hastings and Reginald Rutherford, III 


FLUID AND THERMAL SYSTEMS BRANCH 
PROPULSION DIVISION 

PROPULSION AND VEHICLE ENGINEERING LABORATORY 
RESEARCH AND DEVELOPMENT OPERATIONS 



PRECEDING PAGE BLANK. NOT FILMED 


TABLE OF CONTENTS 

Page 


SUMMARY 1 

4 I. INTRODUCTION 2 

II. BASICS OF LIQUID- SOLID- VAPOR SYSTEMS 4 

III. REVIEW OF PREVIOUS ANALYSES 7 

%> IV. THEORETICAL ANALYSIS ............... 15 

* 

A. Total Energy of a Capillary System ................ 15 

B. Principle of Minimum Energy 20 

C. Application of Variational Principle 22 

V. THEORETICAL RESULTS 29 

A. Bond Number and Container Shape 29 

B . Contact Angle 31 

VI. EXPERIMENTAL VERIFICATION 45 

VII. CONCLUSIONS 51 

APPENDIX A - LIQUID- VAPOR-SOLID INTERFACES 53 

APPENDIX B - BOUNDARY CONDITIONS DEPENDENT ON 
» CONTAINER SHAPE 57 

APPENDIX C - COMPUTER PROGRAM FOR DETERMINING 

INTERFACE SHAPE 64 

REFERENCES 87 


m 



LIST OF ILLUSTRATIONS 


.& 


Figure Title Page 

1 Reference Coordinate System . . , 8 

2 Reference Coordinate System . 16 

3 Reference Polar Coordinate System 1? 

4 Low Gravity Zero Contact Angle Interface Shapes 

in Cylindrical Container s . 33 

5 Zero Contact Angle Interface Deviation from the 

Infinite Bond Number Level in a Cylinder ..... 34 

6 Zero Contact Angle Interface Shapes in Spherical 

Containers for Bond Number =5 35 

7 Zero Contact Angle Interface Shapes in Spherical 

Containers for Bond Number =20 36 

8 Zero Contact Angle Interface Shapes in Shperical 

Containers for Bond Number =70 37 

9 Zero Contact Angle Interface Shapes in Spherical 

Containers for Bond Number = 150 . . * 38 


10 Zero Degree Contact Angle Interface Shapes in an 

Oblate Spheroid 39 

11 Zero Degree Contact Angle Interface Shapes in a 

; Prolate Spheroid . ^0 

12 Effect of Contact Angle on Surface Shapes in a 


Cylinder 41 

13 Contact Angle Influence on Surface Deviation from 

the Infinite Bond Number Position in a Cylinder 42 

14 Effect of Contact Angle on Surface Shapes in a 

Sphere at Bond Number = 0 43 

15 Effect of Contact Angle on Surface Shapes in a 

Sphere at Bond Number =50 . 44 

16 Experimental and Theoretical Interface Shapes 

in a Cylinder . 48 


/ 


* 


IV 



LIST OF ILLUSTRATIONS (Concluded) 


"A 

Figure Title Page 

* 17 Experimental and Theoretical Interface Shapes 

in a Cylinder 49 

* 

18 Experimental and Theoretical Interface Shapes 

in a Cylinder . . 50 

1A Contact Angle Measurement 54 

2A Surface Tensions at a Liquid- Vapor -Solid 

Interface 55 

IB Geometry for Prolate and Oblate Spheroid 

Boundary Conditions 58 

2B Geometry for Boundary Conditions in a Cylinder 62 


% 


V 



DEFINITION OF SYMBOLS 

A Area, ft 2 

a Local acceleration, ft/sec 2 * 

b Horizontal semi -axis of ellipse * 

4 ' 

Bj^ Bond number based on R Q * 

C (y 2 + y' 2 ) 

E Total energy, ft-lbf 

F Integrand of integral to be minimized 

G Integrand of integral to be held constant 

g Q Dimension constant, lb m ft/lb£ - sec 2 

H Bond number based on r Q 

H' Bond number based on r c 

K q Parameter related to curvature at interface centerpoint 

k Q Curvature at interface centerpoint 

L Arc length, ft ^ 

M Mass, lbf - sec 2 /ft 

P. E. Potential energy, ft - lb£ 

P Pressure, lbf/ft 2 

Q Heat, ft - lbf 

R Container radius, ft 

R 0 Characteristic container dimension, ft 

r Radius of curvature, ft 

r — — 

c P _ p- 

■^v o 

S.E. Surface energy, ft-lbf 

S Entropy, ft-lbf/°R 

T Temperature, °R 

U Internal energy, ft-lbf 


vi 



DEFINITION OF SYMBOLS (Concluded) 


ih 


V 
W 
X 

Y 

Yfi 

y' 

y 

y 0 

H 

z 

a 

P 

0 

6i 

\ 

P 


a 

c 

i 

s 

v 

w 


Volume, ft 3 
Work, ft-lb f 
Horizontal distance, ft 

Horizontal distance to container wall at (y, 0) 
Horizontal distance to container wall 
dy/d0 * 

Distance from origin to liquid surface 

Distance from origin to low gravity surface centerpoint 

Distance from origin to surface for Bond Number = ® 

Vertical distance, ft 

Contact angle, degrees 

Empty fraction 

Angle measured from vertical axis to y, degrees 
0 at interface /container wall intersection, degrees 
Lagrange multiplier 
Density, lbf sec 2 /ft 4 
Surface tension, lbf/ft 

SUBSCRIPTS 

Acceleration 

Capillary 

Liquid 

Solid 

Vapor 

Container wall 


vii 



LOW GRAVITY LIQUID- VAPOR INTERFACE 
SHAPES IN AXISYMMETRIC CONTAINERS 
AND A COMPUTER SOLUTION 

By 

Leon J. Hastings and Reginald Rutherford, III 
SUMMARY 

The purpose of this study was to derive a convenient method for 
determining equilibrium liquid- vapor interface shapes in axisymmetric 
containers of arbitrary shape. Several legitimate methods for deter- 
mining interface shapes were previously proposed by other investiga- 
tors. In fact, the basic differential equation for computing interface 
profiles was presented by Bashforth and Adams in 1883. However, 
these solutions are either inconvenient to apply or are restricted to 
a certain range of boundary conditions. 

A convenient form of the basic interface differential equation was 
derived from the familiar principle of minimum surface and potential 
energy using the calculus of variations. The use of a polar coordinate 
system eliminated the convergence difficulties encountered in the previ- 
ous solutions. Also, the derivation enabled the incorporation of a Bond 
number based on a characteristic container dimension into the basic 
differential equation as opposed to a Bond number based on interface 
radius of curvature. 

The basic differential equation and boundary conditions dependent 
on container shape were programmed for a GE235 computer so that sur- 
face shapes for any particular combination of Bond number, vapor 
volume, and contact angle can be determined. The computer solution 
utilizes the Runge-Kutta numerical technique and imposes no significant 
limitations on contact angle or Bond number . 



Representative surface shapes were computed to determine the 
influence of contact angle. Bond number , and container fill level for 
three container shapes: cylindrical, spherical, and spheroidal. It was 
determined that, in a cylinder, the influence of Bond number on inter- 
face deformation is maximum between Bond numbers of two and twenty 
and becomes negligible for Bond numbers greater than approximately 
200. In spherical or spheroidal containers, the empty fraction has a 
significant effect on the interface profile. Also, unlike the cylinder, 
a contact angle of 90 degrees does not assure negligible interface 
distortion in spherical or spheroidal containers. The limiting contact 
angle is that angle measured in the liquid between a horizontal plane 
corresponding to the infinite Bond number liquid level and the tangent 
to the container boundary. 

The theoretical profiles were compared with experimental data, 
and exceptional agreement was obtained. In fact, if the actual contact 
angle is known, it is believed that the static equilibrium interface pro- 
files can be computed with greater accuracy than they can be measured 
due to the distortion and reflection problems inherent in such experi- 
mental measurements. . 

I. INTRODUCTION 

In environments devoid of any disturbances except that of a low 
acceleration or gravity, surface tension forces become comparable to 
those of gravity, and equilibrium liquid-vapor interface shapes may 
radically depart from the near flatness observed in normal gravity. 
This is especislly true of the many liquids that exhibit wall contact 
angles at or near zero degrees. The problem of low gravity interface 
descriptions for various container shapes, liquids and acceleration 
levels is of especial interest to engineers responsible for the design of 
propellant control schemes for space vehicles and storage tankers that 
must operate for long periods of time in orbital environments, that is, 
under low gravity conditions. 



The purpose of this study was to provide a general method that 
can be used to determine low gravity interface shapes for liquids in 
axisymmetric containers subjected to axial accelerations less than that 
of normal gravity. As demonstrated herein, a single differential equa- 
tion, which is applicable to all containers that are symmetric about an 
axis parallel to the acceleration direction, can be derived from the 
principle of minimum potential and surface energy using the calculus of 
variations. This differential equation in conjunction with boundary con- 
ditions dependent on container shape and contact angle can be solved 
numerically for the desired interface shape using the Runge-Kutta 
iteration technique. Although the basic differential equation and its solu- 
tion is applicable to all containers symmetric about the vertical axis, the 
following geometrical shapes were selected for analysis based on their 
practical significance in space vehicle applications: (a) cylindrical, 

(b) spherical, and (c) prolate and oblate spheroids. This theoretical 
method for interface determination was programmed for a GE 235 com- 
puter, and theoretical interface shapes were computed for a wide range 
of conditions. .Since practically all known liquid propellants considered 
for space vehicle propulsion exhibit contact angles at or near zero 
degrees, interface configurations for fluids of zero contact angle are 
emphasized. 

Also presented are experimental data concerning low gravity 
interface shapes that were obtained from a Lockheed Missiles and 
Space Company experimental program. Although the data are not 
extensive in scope, it is sufficient to substantiate the proposed theoreti- 
cal techniques of interface shape determination. 
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II. BASICS OF LIQUID- SOLID- VAPOR SYSTEMS 


A fundamental property of liquid surfaces is their tendency to 
contract to the smallest possible surface area for a given volume, i. e. , 
a spherical surface. This property has been proven by thorough 
theoretical and experimental investigations which have been performed 
ever since the phenomena of capillarity was first noted by Leonardo da 
Vinci (Ref. 1)*. The properties of molecules in liquids can easily be 
used to account for this minimization of surface area. On the interior 
of the liquid, each molecule is surrounded by others on every side and 
is, therefore, subject to attraction in all directions. This, however, 
is not the case in the surface. Molecules in the surface are attracted 
inwards and to each side by neighboring molecules, but there is no out- 
ward attraction to balance the inward pull since the molecular density 
of the vapor is much less than that of the liquid. Hence, every surface 
molecule is subject to a strong inward attraction perpendicular to the 
surface. This inward attraction causes the surface molecules to move 
inwards more rapidly than others move outwards to replace them. The 
diminishing of molecules in the surface continues until the maximum 
number of molecules are in the interior, that is, until the surface is 
spherically shaped, subject to external conditions or forces acting on 
the liquid . 

When external forces, such as gravity and those created by the 
presence of other materials, must be considered, the resulting liquid- 
vapor interface shapes are more complicated. Gravity forces always 
tend to minimize the potential energy associated with the mass of 


* This reference contains an excellent summary of early works per- 
formed on capillary phenomena. 
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liquid affected by the liquid- vapor interface shape and, therefore, is an 
important factor in determining equilibrium interface shapes. In addi- 
tion, another influencing factor that must be considered is the presence 
of a solid material and the molecular attraction between it and the 
liquid. If the molecules in the liquid are attracted by the solid (adhesive 
forces) more than by neighboring molecules in the liquid (cohesive 
forces), the liquid is said to "wet" the solid, that is, the liquid tends to 
spread on the solid surface. If the cohesive forces are greater than 
the adhesive forces, the liquid does not wet the solid surface. 

Since surface energy is associated with a liquid- vapor interface 
and work must be done against the internal liquid molecules to extend 
the surface, a vast number of problems relating to the equilibrium 
position of surfaces can be solved if the magnitude of surface energy is 
known. To simplify the calculations, however, a hypothetical tension, 
which acts in all directions parallel to the surface, is substituted for the 
surface energy. This hypothetical tension is generally termed "surface 
tension." Surface tension has the same dimensions as surface energy 
per unit surface area and it must have the same numerical magnitude. 
Proof of this is relatively simple and is readily available in literature 
on physics of surfaces (a thermodynamic proof is contained in Ref. 2). 
Several references (see References 3 and 4) emphasize that the concept 
of liquid surfaces behaving like a stretched membrane must not be 
misconstrued, since surface energy is the fundamental liquid pro- 
perty and surface tension is merely its mathematical equivalent. 

At any rate, the proper utilization of surface tension is very con- 
venient. In systems involving liquid surfaces the equilibrium position 
can be acquired by totaling up the changes in surface energy of the 
various interfaces (liquid- solid, liquid-vapor, and solid- vapor) whose 
areas are altered by a displacement. However, if the surfaces are 
considered at the boundaries where each is pulling with the appropriate 



"surface tension" on their boundaries, then the changes in liquid-solid 
and solid- vapor surface areas do not have to be calculated. Such an 
approach is taken in this report, which considers the effects of low 
gravity environments, container shapes, liquid "wettability"*, and 
liquid surface tension, on the equilibrium interface shape. 


* The term "wettability" refers to the degree of solid-liquid attrac- 
tion. The degree of wetting is usually measured in terms of "con- 
tact angle, " which is explained in Appendix A. 




in. REVIEW OF PREVIOUS ANALYSES 


Several investigators have devoted considerable time to the 
evaluation of static equilibrium configurations of liquid-vapor systems 
in zero and low gravity environments. Most of their evaluations, how- 
ever, apply only to containers of cylindrical shape (see for example 
References 5, 6, and 7). Only the works of Bashforth and Adams 
(Ref. 8); Reynolds, Saad, and Satterlee (Ref. 9); and Li (Ref. 10) can 
be applied to axisymmetric containers of arbitrary shape. 

Undoubtedly, the most important early exploration of the subject 
was performed by Bashforth and Adams and is recorded in their book 
which was published in 1883. Bashforth and Adams were primarily 
concerned with the shape of liquid droplets on or suspended from a 
horizontal surface. Their results apply reasonably well to the 
present problem as demonstrated by Yeh and Hutton (Ref. 11) and 
Jurney (Ref. 12). 

Bashforth and Adams derived the governing differential equation 
beginning with the generally accepted relation for pressure difference 
across any curved liquid-vapor interface at a particular point 

p v - P * ■ °iv + 77 ) < l 

where P v and Pjg are the pressures on the vapor and liquid side of the 
surface respectively, r I and r 2 are the principal radii of curvature of 
the surface at the point of interest, and (rg v is the surface tension of 
the interface. 

As shown in Figure 1, if Pjp 0 is the interior pressure of the sur- 
face at its origin, then due to hydrostatic pressure variation 
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Note : Tj is in the XZ plane and 

r 2 is in a plane perpendicular 
to the XZ plane. 



FIGURE 1. REFERENCE COORDINATE SYSTEM 
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p i = p f o - Pi a z 


( 2 ) 


where is the liquid density and a is the local acceleration. 
Also, at the origin 


p _ P _ t-JLi.Y. 


(3) 


where r j = r 2 = r Q in equation (1). 

Using equation (3) and substituting for P from, equation (2), 


_ 2 erg - 


p v - p i * “ + 


Z = <r lv +■—) 


(4) 


Let X be the horizontal and Z be the vertical coordinates of any 
point in a meridional section of the surface, rj the radius of 
curvature of the meridional section at that point, and <|> the angle which 
the normal to the surface makes with the axis of revolution. Then, the 


length of the normal terminated by the axis (r 2 ) is 
(4) becomes 


X 


sin cj> 


, and equation 


^ + ^SJfe = 2 + H (*) 
*i X/r 0 Vr 0 / 


(5) 


where H = - r<? is a dimensionless number often termed the Bond 

°iv 

number . 

Another form of equation (5) can be easily derived. Since the 

expression for the radius of curvature of a line in the XZ plane is 

2 1 s/z 


r, = 


1 + 


(— V 

\dX/ 


d z Z 

dX 2 


and from Figure 1 
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then equation (5) is equivalent to 


d 2 Z 
dX 2 + 



dZ / 2 HZ\ f /dZ\ 2 1 3 / 2 
XdX = \r Q + r 0 2 / [_ 1 + \dX/ 


(6a) 


or 


_L_ _d_ 
X dX 


X 


1 + 


dZ 

dX 

(— y 

VdX/ 


172 


> * 


_ 2 _ 

r„ 


HZ 


(6b) 


These equations will be recognized in subsequent paragraphs as forms 
of the differential equation derived by other investigators using other 
approaches to the problem. 

Equation (5) was solved by Bashforth and Adams using either the 
arc length L or <j> as the independent variable. For example, if 
cj> is taken as the independent variable, then upon integration of the 
equations 

dX . dZ . , 

— = rj cos 4>, — = r x sm<j> (7) 

with use of the initial conditions at the origin 
«j> = Z = 0; rj = r 2 = r Q 

equation (5) will yield r j as a function of the coordinates X and Z. 

First, the form of the curve in the neighborhood of the origin can 
be found by developing r x and the coordinates X and Z in series of 
ascending powers of cj). Next, the coordinates for larger values of <j> 



can be obtained by step by step numerical integrations. Extensive 
numerical tables are presented in Reference 8 that give coordinates 
of surfaces for various values of H, 

In present day space vehicle applications, it is desirable to base 
the interface shape on a Bond number in which a container dimension, 
such as tank radius, is the characteristic dimension. Bashforth and 
Adams' system can be converted to the desired Bond number system, 
but extensive interpolations are often required between tables in order 
to locate the proper liquid volume, contact angle (a), and Bond number 
combination. In addition, according to Yeh and Hutton's studies, 
Bashforth and Adam's tables contain a maximum Bond number (based 
on container radius) of about ten, so that for many applications their 
tables must be considerably extended. 

Reynolds, et al obtained a modified but equivalent form of 
Bashforth and Adams' differential equation by performing a force 
balance on an infinitesimal annular ring cut from the meniscus. Using 
the coordinate system shown in Figure 1, a vertical force balance on 
an annular ring yields 

Capillary Force, Fc = Pressure Force, Fp 
or 

d (2 ttX o -i v ||) = (P v - P i ) (2-rr XdX) (8 ) 

Substituting equation (2) for P^ , 
d / dZ\ XdX 

°iv dL \ dL/ dL ^ Pv P fo + P a ) (9) 

Introducing the dimensionless quantities 
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H i = £L*L sl 
^iv 


where r c = 


JUjL. 


(Py - Pi o> 


the final differential equation becomes 


d / dz\ dx TTt . 

S V x di) = x dI ll + H ' z) 


( 10 ) 


A second differential equation was obtained from the geometric condition 


(If) + df) = 1 


which when differentiated, gives 


dz d 2 z dx d 2 x _ 

df df 2 + df dT 5 " 


(ID 


Equations (10) and (11) form a pair that can be solved numerically. 
The problem is treated as an initial value problem in which 


x (0) = z (0) = z'(0) = 0 

and 

x' (0) = 1 

Calculations of this type were performed for various values of H' 
(from 0 to 10) and contact angles. The calculations are similar to those 
of Bashforth and Adams and are subject to similar limitations. 

In fact, it can be shown that equation (11) is a modified form of 
Bashforth and Adams’ equation. Beginning with equation (9) and recall- 
ing that 



r 
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Pjg aZ HZ 


and. 


-=M§)T ax 


then equation (9) becomes 

dZ 

I X 

J_ _d_ 

X dX 


“ /dZ\ 2 1 1 / 2 

1 + (dx) 


_2_ + HZ 


which is identical to equation (6b). 

Using a still different approach, that is, the principle of minimum 
surface and potential energies and the calculus of variations, Ta Li* 
derived the interface differential equation 



where X is a Lagrange multiplier; and atZ = 0 


d 2 Z _ J_ dZ 
dX* = r Q * dX 


0 , and Z s Zq " 0 


Then Li's final equation can be written 

,-i / 2 | 


X 


d 

dX 


X 


dZ 

dX 


/ dZ \ 2 ' ' 

1 + (dx) _ 


HZ 
= „ 2 


which is identical to equation (6b). 


* Li's derivation is not presented here since a similar procedure is 
used in the subsequent section of this report. 
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However, Li solved the differential equation through a transfor- 
mation of variables and the introduction of a power series to give 

“ = f < X > 
r o 

Li obtained the first four terms of the expansion, and later Yeh and 
Hutton derived and used three additional terms of the expansion. How- 
ever, the minimum contact angle for which Yeh and Hutton obtained 
satisfactory convergence was about 53 degrees. Yeh and Hutton con- 
cluded that Li's solution could not be used for contact angles at or near 
zero degrees, since the number of coefficients required for the series 
solution was a function of the cotangent of the contact angle, and there- 
fore, an infinite number of coefficients was required as the contact 
angle approached zero. 

In conclusion, it is notqd that several approaches can be used to 
derive the basic differential equation for low gravity surface profiles. 
However, the independent variables chosen, the method of solution, and 
the coordinate system are important factors in obtaining numerical 
results for a wide range of conditions. 



IV. THEORETICAL ANALYSIS 


As related previously, Ta Li has proposed a general method of 
calculating low gravity liquid- vapor interface shapes in axisymmetric 
containers. His method appears to be legitimate but is unsuccessful in 
the many cases when the interface slope reaches or passes the vertical, 
i.e. , curves back on itself. Tq avoid this difficulty, a polar coordinate 
system is proposed in the following sections. An integral equation is 
derived from energy considerations, and the calculus of variations is 
used to obtain a differential equation that can be solved numerically. 

A. Total Energy of a Capillary System 

Consider a container that is cylindrically symmetric and is partly 
filled with liquid. The total energy of this arrangement is the sum of 
surface and potential energies. 

Total Energy = Surface Energy + Potential Energy 

E — - S . E . 4* P . E . 

The surface energy is the sum of surface energy of three interfaces: 
liquid-vapor (fv), solid-vapor (sv), and solid-liquid (si)* Therefore, 

S. E. = v Ajg v + c sv A gv + <r S £ A g ^ 

where <r is the surface tension and A is the area of the respective inter- 
faces. Using the fact that the total area of the solid surface A^q^j is 
constant and the relation* 


o-j[ v cos a s or sv - cr s ^ 

where a is the wall/liquid- vapor interface contact angle, then 


* Refer to the Appendix A for a discussion of this relation. 
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S.E. — crj£ v (Ajp v - A Sj g cos o;) + o"s v ^total 


(12a) 


or 


S. E. = <r£ v (Aj{ v - A Sj j cos a) + Constant 


(12b) 


Consider the coordinate system, indicated in Figure 2. The independent 
variable is Z, the liquid- vapor surface is described by X = X(Z); and 
the container boundary, which is axi s ymme tr i c about the Z axis, is 
described by X w = X W (Z). 

Z 



Acceleration 

Direction 


X W (Z) 


FIGURE 2. REFERENCE COORDINATE SYSTEM 
The liquid- vapor interface surface area is 

dA iv = 2 t rX (dZ 2 + dX 2 ) l/2 (13) 

and the liquid volume (V£), which is constant, is 

dV£ =tt(X w z - X 2 ) dZ (14) 
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Also, the potential energy of the liquid referenced to an arbitrary point 
on the liquid surface is 

d(P. E. ) = a( Pj? - p v ) Z dV£ 
or 

d(P. E. ) =traAp (X w 2 - X 2 ) ZdZ (15) 

where: a = local acceleration of gravity 

Ap = liquid density - vapor density = - p v (usually 

Ap = p^ for all practical purposes) 

Now, consider a change to the coordinate system indicated in Figure 3. 



FIGURE 3. REFERENCE POLAR COORDINATE SYSTEM 

The parameter y describes the liquid-vapor interface shape, and the 
angle 0, measured from the vertical to y, is the independent variable. 
The coordinate transformation is 


* 

Z = - R 0 y cos 0 

(16a) 


X = R c y sin 0 

(16b) 


where R 0 is any convenient length parameter of the container. The 


17 



container boundary is now given by 


Y = f (y, 0) 


(17) 


where Y is the dimensionless radius of the container wall in a horizontal 
plane with the point (y, 0). 

In the new coordinate system, if relations (16) and (17) are 
substituted into equations (13), (1,4), and (15) 

T /dy\ 2 ”]*^ 2 

d Af v = 2ir R 0 2 y sin0 y 2 + J d0 (18) 


d Y a * ttR 0 3 (Y 2 - y 2 sin 2 0)(y sin 0 - dy/d0 cos0)d0 (19) 

d(P„E.) = -ir aApR 0 4 y cos 0 (Y 2 - y 2 sin 2 0)(y sin0 


- dy/d0 cos 0)d0 (20) 


Note that the total liquid volume must remain constant and from equa- 
tion (19) is 


0 i 

Vj? =irR 0 3 ^ (Y 2 - y 2 sin 2 0)(y sin0 
0 


- y' cos0)d0 + f (Z Q ) 


( 21 ) 


where f(Z Q ) is the volume of liquid below the horizontal plane at 0 a 0 
and is, therefore, independent of the detailed form of interface shape. 
Using relations (12) and (18), anewrelation can be derived for surface 
energy. Referring to Figure 3 and letting 0j be the final 0, i.e., 0 
measured at the interface /container wall intersection, 

9i 


y — 2 IT R-£ 


I 


ysin0(y z +y ,2 ) 1/z d0 


\i/z, 


( 22 ) 


0 


Note: y' 


dy 

d0 
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Substituting equation (22) into (12) 

0i 

S.E. = 2irR 0 2 ar lv £ J y sin 0 (y 2 + y ' 2 ) d0 - — ^ ^ 2 ~ J 


+ Constant 

However, A SJ j is independent of the detailed form of y(0), the interface 
shape, and is dependent only on 0 1 and the container boundary. 
Therefore, 

-nA = «<».> 

TT K.Q 


and 


0i 


S.E. = 2Tr(r iv R 0 2 | J y sin 0 (y z + y' 2 ) l/z d0-f(0 1 )^osa 

0 

+ Constant 

From equation (20), the total potential energy is represented by 

0i 

P. E. = - tt a A p R 0 4 y y cos 0 (Y 2 - y 2 sin 2 0)(y sin0 - y* cos 0)d0 


(23) 


(24) 


If a dimensionless parameter termed "Bond number," which is a ratio 
of "body forces" to "surface tension forces," 

ApaRo 2 

b N = 


o'fv 


is incorporated, then equation (24) becomes 

0i 

i 

P. E. = - ^ 

0 


j^BnRq 2 y y cos 0 (Y 2 - y 2 sin 2 0)(y sin0 

- y' cos 0) d0 


(25) 
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By combining equations (23) and (25), the total energy is represented by 

E = iro‘ J j v R 0 2 J^2 y sin0 (y 2 + y' 2 ) 1 / 2 d0 - 2f(0j) cos a 

8 , ° < 26 > 

- B N J y cos 0 (Y 2 - y 2 sin 2 0)(y sin 0 - y' cos 0)d0^ + Constant 


In subsequent sections the expressions derived for total liquid 
volume, equation (21), and total energy, equation (26), will be used in 
conjunction with certain physical considerations to provide a basic 
equation for the equilibrium liquid- vapor interface profile. The physi- 
cal considerations are that the total liquid volume must remain constant 
and that the principle of minimum energy must be satisfied. 


B. Principle of Minimum Energy 

Numerous investigators have observed through experimentation 
that the total energy, E, of a liquid-vapor- solid system such as 
the one under consideration, tends toward a minimum. The thermo- 
dynamics of such a capillary system can be used to explain this 
principle of minimum energy. Using equation (12b) the work required 
to extend the liquid-vapor interface is given by 

d W c = d(S. E, ) = <Tjg v dA c 
where A c = A| v - Af s cos a. 

The work performed against gravity when the capillary system's 
center of mass is moved upward a distance of dZ, is 

d W a * MadZ . 

Then the total reversible work done on a capillary system is 

d W = o-jg v dA c + MadZ (27) 
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In a reversible or irreversible change of state, the entropy 
change, dS, is related to the energy transferred to the capillary sys- 
tem by 

T dS = dU - dW (28) 

where U is the internal energy and T is the temperature. 

Substituting equation (27) in equation (28) 


T dS = dU - q-jg ydA c - MadZ (29) 

The general criteria for neutral equilibrium of a system for all 
possible variations is 

AS)g = AE)g = 0 

and for stable equilibrium is 

AS) g < 0 or AE)g > 0 

where £ is the stored system energy. Thus, in general, for all possible 
variations in the system for which 


AS) e £ 0 or A E) s £ 0 

a stable or neutral equilibrium state exists. 

In this case, equation (29) represents the infinitesimal change 
between two equilibrium states. Consider changes at constant internal 
energy ( £ = U in this analysis) 


dS)u = - 


. V £ v^^c 
T 


MadZ 

T 


(30) 


Then, if the initial state is in neutral or stable equilibrium and 
since T is always positive 


dS)jj | 0 or cr^ v dA c + MadZ 2 0 

and 

dE > 0 
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Since the capillary system energy can only increase, the equilibrium 
state must be one of minimum energy , E. 

Thus, the energy described by equation (26) must yield a mini- 
mum in order to satisfy the principle of minimum energy. As illu- 
strated in the next section, this principle of minimum energy enables 
a convenient application of the calculus of variations . 

C. Application of Variational Principle 

In view of the aforementioned factors concerning the total energy 
and liquid volume of the system, it is necessary to minimize the total 
energy while keeping the liquid volume constant. Such a condition is 
known as an "isoperimetric problem" with a mobile upper limit (0j) 
in the calculus of variations. As shown by Li, it is convenient to use 
the calculus of variations to construct a function from Euler's relation 
(Ref. 13) 


9(F + XG) 

8y 


_d_ r a(F + xg) ~ 
de [ V 


(31) 


where "F" and "G" are integrands of integrals to be minimized and 
held constant, respectively, and "X." is a Lagrange multiplier. 
Ther efore, 

F = F (y, 0)* 

= 2 y sin 0(y 2 + y' 2 ) l ^ z 

- Bjvj y cos0(Y 2 - y 2 sin 2 0)(y sin0 - y' cos 0) 


and 

G = G (y , 0) 

* The term £{0^ is not included as a part of F because, as proven by 
Li, this only effect of the term is to assure the proper contact 
angle at the container boundary. Thus, this terminal condition 
is not satisfied for the time being. 
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G = (Y 2 - y 2 sin 2 0)(y sin0 - y' cos 0) 


Performing the indicated operations in equation (31), 

“ (F + XG) = 2 sin 0 (y 2 + y ,2 ) l/z + 2y 2 sin0(y 2 +y' 2 )- x ! z 

- cos 0(Y 2 - y 2 sin 2 0)(y sin0 - y' cos 0) 
r9Y 2 


foY \ 

\~ Qy ~ - 2 y sin 2 QJ (X - Bfj y cos 0) (y sin 0 - y' cos 0) 


+ sin0 (X - Bj^ y cos0)(Y 2 - y 2 sin 2 0) 


(32) 


and 


d 

d0 


8(F + XG)" 

_d_ ' 

9y * 

= de 


2yy' sin0(y 2 +y ,2 ) -1 / 2 J 


+ (Bjvjy cos 0 - X)(Y 2 - y 2 sin 2 0) cos 0 

e 2y y n sin0(y 2 + y ,2 ) _1 / 2 +2y ,z sin0(y 2 +y ,2 ) _1 / 2 
- 2yy' sin0(yy' +y'y")(y 2 + y’ 2 ) -3 / 2 

+ 2y y 1 cos 0(y 2 + y ,Z ) ” 1 / 2 - (B]vj y cos 0 - X)(Y 2 - y 2 sin 2 0) sin 0 


+ (Bj^ycos 0 - X) 


-jt (Y 2 ) - 2y y ' sin 2 0 - 2y 2 sin 0 cos 0 Icos 0 
d0 J ' 


) < 


+ Bj«j (y* cos 0 - y sin 0)(Y 2 - y 2 sin 2 0)cos 0 


(33) 


On combining relations (32) and (33) according to equation (31) and with 
considerable algebraic manipulation, it can be shown that: 

0 = (2y 2 +3y' 2 -yy.')C-’/ 2 - - C ~ ^ S?* 8 

+(X.B N ycos8)[^(^-+g^)-l] (34) 


23 



where C = (y 2 + y' 2 ) for simplification 


On examining 


8Y 2 cote 9Y 2 
9y y 9 0 


where Y =f(Z)and 


Z a f (y, 0) a - R q y cos 0 


Note that 

8Y 2 dY 2 9Z 
8 y dZ 9y 

9Y 2 dY 2 8Z 
9 0 " dZ 9 0 


a - Rq COS 0 


dY 2 

dZ 


= + R 0 y sin 0 


dY 2 

dZ 


Hence, 

9Y 2 cot 0 9Y 2 

9y y 90 


(35) 


This relation is very significant, for it means that container shape 
drops out of the differential equation (34). Incorporating equation (35) 
in equation (34) and solving for y" 

n 2 y 2 + 3 y' 2 y 1 cot 0 (y 2 + y* 2 ) 


+ ^ (Bj^ y cos 0 - X.)(y 2 + y > 2 ) 3 / 2 


(36) 


The imdeter mined multiplier X. can be solved for by treating the 
problem as an initial value problem in which the following "initial 
conditions" are prescribed 

y(o) =Yo 

y' (0) = 0 (interface symmetry) 
y" (0 ) =y Q (1*- K 0 ) 

The initial value for y"(0) can be derived using the equation for 
curvature of an arc at any particular point 
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Substituting the conditions specified for y(0) and y'(0) and solving for 
y"(0) 


Y"(0) = y Q ( - y 0 k 0 ) 

or 


y"(0) =y 0 d - k q 

where K 0 = y Q k o and may be considered a parameter related to inter- 
face curvature at the point of symmetry. Therefore, by substituting 
the prescribed "initial conditions" in equation (36) and solving for X 

x 2 K 0 „ 

X = ~y^ + Bn Yo 

Using this value for X, the final form of the differential equation is 
obtained 

y u - if. ± 3y! , i . xl. cote(y 2 +y .2) 

+- f ® N (y c°se - y 0 > - -^lly 2 + y ,2 ) 5/2 (37) 

y l ■ yo j 

The validity of equation (37) can be checked at the two limiting condi- 
tions, namely, = 0 and as Bj^ approaches infinity. At B-^ = 0 (zero 
gravity), it has been experimentally verified by numerous investigators* 
that the surface of a totally wetting liquid will tend toward a shape of 
minimum surface area, that is, a sphere. Therefore, at Bjj = 0 equa- 
tion (37) must be satisfied by the equation for a circle, y = y Q cos 0. 


* See, for example. Reference 14. 

25 



Substitution of 


y = y 0 cos 6 
y' = - y 0 sin0 

K 0 = l-^ 

° y Q 


+ 2 


in equation (37) gives 


y" = - y Q cos 0 

which verifies that y = y Q cos 0 at Bjsj = 0 . 

As Bjsj approaches infinity the equation should become that for a 
horizontal line. That this occurs is easily verified by dividing equa- 
tion (37) by Bjsj and taking the limit as Bj^ approaches inifinity. In such 
a case, equation (37) becomes 


0 = 


_1 

y 


(y cos 0 - y Q ) 


or 

y Q 

^ - cos 0 


which is the equation for a horizontal line. 

For values of Bjsj between the two end points, numerical solution 
of equation (37) is easily accomplished using the Runge-Kutta iteration 
technique (Ref. 15), which yields a dimensionless plot of the liquid- 
vapor interface for a given y Q , K 0 , and Bj^. However, the additional 
boundary conditions of contact angle, constant vapor (or liquid) volume, 
and container shape must be satisfied in order to obtain the desired 
interface shape. The equations for calculating these boundary condi- 
tions are all functions of container geometry and, therefore, constitute 
the only changes that must be considered when the container shape is 
modified. To illustrate the modifications that are required by a change 
in tank configuration, the equations for boundary conditions dependent 
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on tank shape are described in Appendix B for four shapes: spherical, 
prolate spheroid, oblate spheroid, and cylindrical. 

Equation (37) and the boundary conditions described in Appendix 
B are programmed for a GE 23 5 computer so that the surface shape 
for any particular combination of Bond number, vapor volume, contact 
angle, and container shape can be determined. Basically, the calcula- 
tion of a surface shape using the computer program consists of initiat- 
ing the computation procedure with an initial set of conditions at the 
interface centerpoint and using the Runge-Kutta procedure until the 
calculated surface intersects the container wall. Whenever the surface 
reaches the container boundary, the boundary conditions of contact 
angle and vapor volume are checked. If these boundary conditions are 
not satisfied, appropriate changes in the initial set of conditions are 
selected by the computer and used to initiate the Runge-Kutta procedure 
again. Thus, the iteration procedure continues until all conditions are 
satisfied and the desired interface is yielded. 

The computer program is described briefly in engineering terms 
in Appendix C. This information provides a better visualization of the 
relationships of the various parameters involved in the determination 
of an interface shape. 

Thus, a differential equation (equation (37)) which enables the 
elimination of difficulties incurred in the application of previously 
derived interface differential equations (see Chapter III) has been 
formulated. Desirable features that have been incorporated in equation 
(37) include: 

1. The Bond number, Bj^, appearing in the differential equation 
is based on a characteristic container dimension rather than on the 
surface radius of curvature. 

2. The polar coordinate system utilized eliminates the possi- 
bility of "double -valued" functions which can occur in other coordinate 
systems whenever the surface profile curves back on itself, that is. 



whenever two values of Z occur for a given value of X. Thus, solutions 
are possible for zero contact angles and zero Bond numbers. 

3. A computer solution of the differential equation is easily 
accomplished using the Runge-Kutta numerical technique. Bond num- 
ber, contact angle, container shape, and fill level (except in the case 
of a cylinder) are the only necessary input. 

Theoretical interface shapes for various conditions were deter- 
mined using the computer program. The results are presented and 
discussed in the subsequent chapter. 



V. THEORETICAL RESULTS 


Low gravity interface shapes in spherical, prolate and oblate 
spheroidal, and cylindrical containers for a wide range of Bond num- 
bers, fill.levels, and contact angles were determined and are presented 
in Figures 4 through 15. However, since almost all known liquid pro- 
pellants for space vehicle applications have contact angles near zero 
degrees, the presentation of interface shapes emphasizes shapes for 
zero degree contact angle fluids. A discussion of the influences of 
various parameters is contained in subsequent paragraphs. 

A. Bond Number and Container Shapes 


Figures 4 through 11 illustrate the effects of Bond number and 
tank configuration on surface shapes for zero degree contact angle 
fluids. As one would anticipate, the liquid-vapor interface shape 
approaches the shape observed in normal gravity as the Bond number 
increases. Some specific observations on the interface shapes calcu- 
lated for each of the four container shapes analyzed are as follows: 

1. Cylindrical Containers - Since container empty fraction has no 
influence on interface shape in a cylinder (provided the container top or 
bottom does not interfere with interface formation). Bond number is 
sufficient to prescribe interface shape for a given contact angle. Bond 
number as defined herein for a cylinder is 
,2 


b n = 


PaR 

o'fv 


where the container radius, R, is the characteristic dimension. The 
vertical and horizontal coordinates of the interface shapes presented 
are, therefore, nondimensionalized with respect to container radius as 
follows 


_Z _ Vertical distance measured from Bond number = qq position 
R Container radius 
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X Radial distance from container axis 

R “ Container radius 

Interface shapes for specific Bond numbers ranging from 0 to 200 
are illustrated in Figure 4. However, the influence of Bond number on 
interface shape is best illustrated in Figure 5 where the vertical posi- 
tion of the low gravity interface at the two end points (at the container 
center and wall) are presented as a function of the Bond number from 
one to 1000. For Bond numbers less than two and greater than twenty, 
the interface shape begins to asymptotically approach maximum and 
zero deviation from a horizontal interface. Since a point of inflection 
occurs between Bond numbers of two and twenty, the Bond number 
influence on interface deformation is maximum in this region. Also, 
it is interesting to note that the interface distortion becomes very 
small or even negligible for Bond numbers greater than approximately 
200 . 

If the interface shape is desired for a Bond number that is not 
presented. Figure 5 can be used to determine the position of the inter- 
face at its two end points. These two dimensions will aid in interpo- 
lating the interface data presented to find the required interface shape. 

2. Spherical Containers - Since spherical containers are often con- 
sidered for space vehicle applications, a rather extensive presentation 
of interface shapes is contained in Figures 6 through 9 for Bond numbers 
ranging from 5 to 150. The interface shape for the limiting case of 
Bond number equal zero was omitted because the interface merely 
assumes the shape of a sphere with a volume equal to the vapor volume 
in the container. The dimensionless parameters and symbols used in 
the figures presented are the same as those described for a cylinder. 

The empty fraction has a very significant effect on the interface 
shape in a sphere for a given Bond number. This observation is 
evidenced by comparing the low gravity interface shapes for normal 
gravity liquid levels in the upper half of the container with those in the 
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lower half. The interface shapes in the lower half have significantly 
less curvature than those in the upper half because the interface must 
bend less to become tangent to the container wall and satisfy the zero 
degree contact angle condition. 

3. Oblate and Prolate Spheroids - In cryogenic space vehicle applica- 
tions, particularly liquid hydrogen (fuel)/liquid oxygen (oxidizer) sys- 
tems, the oxidizer tank is usually an oblate spheroid or some modifica- 
tion thereof. As a typical spheroid, the oxidizer tank shape used on the 
Centaur space vehicle was chosen for analysis (1 by 1.38 ellipse). 
Interface shapes for four different fill levels and Bond numbers of 5, 

20, 50, and 100 are shown in Figures 10 and 11, The dimensionless 
parameters used are the same as those for a sphere except that the 
characteristic length parameter, R 0 , is one-half the vertical height of 
the spheroid. The selection of vertical height instead of container width 
as the characteristic dimension was necessary to simplify the integra- 
tion of spheroidal shapes in the computer program. However, whether 
width or length is used in the Bond number is somewhat arbitrary as 
long as care is taken to maintain consistency when the influence of 
Bond number is discussed. 

Examination of the interface shapes presented discloses that the 
interface characteristics are a combination of those noted for cylindri- 
cal and spherical containers. As one would probably anticipate, in 
addition to the effects of Bond number, contact angle, and empty frac- 
tion, the ratio of major axis to minor axis must be considered a very 
significant parameter when determining interface shapes in spheroidal 
container s . 

B. Contact Angle 

As mentioned previously, most liquid propellants considered for 
space applications seem to exhibit zero or near zero degree contact 
angles on solid materials. However, to illustrate the effect of contact 
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angle, representative interface shapes for contact angles from. 5 to 90 
degrees are presented in Figures 12 through 15 for cylindrical and 
spherical containers. 

1. Cylindrical Containers - Figure 12 illustrates the effect of contact 
angle on interface profile for Bond numbers of 0 and 50; and as one 
would expect, the influence of contact angle decreases with increasing 
Bond number. The influence of contact angle is best demonstrated in 
Figure 13, where the interface rise above the infinite Bond number 
liquid level versus Bond number for various contact angles is presented. 
For example, examination of this figure reveals that the difference 
between surface shapes with 0 and 5 degree contact angles becomes 
almost negligible near a Bond number of 100. 

2. Spherical Containers - High contact angle surface profiles are pre- 
sented for three fill levels and Bond numbers of 0 and 50 in Figures 14 
and 15. Unlike the cylinder, a contact angle of 90 degrees does not 
assure negligible interface distortion at all Bond numbers. The limit- 
ing contact angle in all container shapes is that angle measured in the 
liquid between a horizontal plane corresponding to the infinite Bond 
number liquid level and the tangent to the container boundary. Hence, 
in all vessels with curved boundaries, this limiting angle is dependent 
on fill level. It is not surprising, therefore, that for the fill level near 
the bottom of the sphere, the interface becomes flat for all Bond num- 
bers when the contact angle is equal to 54 degrees. At the fill level 
near the top of the sphere a contact angle of about 143 degrees is 
required before Bond number no longer affects the surface shape. 
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FIGURE 4. LOW GRAVITY ZERO CONTACT ANGLE INTERFACE 
SHAPES IN CYLINDRICAL CONTAINERS 
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FIGURE 6. ZERO CONTACT ANGLE INTERFACE SHAPES IN 
SPHERICAL CONTAINERS FOR BOND NUMBER = 5 
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FIGURE 8. ZERO CONTACT ANGLE INTERFACE SHAPES IN 

SPHERICAL CONTAINERS FOR BOND NUMBER = 70 
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FIGURE 9. ZERO CONTACT ANGLE INTERFACE SHAPES IN 

SPHERICAL CONTAINERS FOR BOND NUMBER = 150 
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FIGURE 12. EFFECT OF CONTACT ANGLE ON SURFACE 
SHAPES IN A CYLINDER 
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FIGURE 14. EFFECT OF CONTACT ANGLE ON SURFACE 
SHAPES IN A SPHERE AT BOND NUMBER = 0 
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VI. EXPERIMENTAL VERIFICATION 


Experimental surface shapes for various Bond numbers and 
liquids in cylindrical containers were measured in normal gravity 
by personnel at Lockheed Missiles and Space Company, and are 
described in Reference 16. The Lockheed Company data are used 
herein to verify the theoretical solution presented in this report. 

The test equipment and procedure used by Lockheed was 
relatively simple. The test containers were fabricated from Lucite 
blocks approximately 2x2x4 inches by drilling and polishing holes of 
various diameters in the blocks. The experimental procedure consisted 
simply of placing the test liquid in containers of various diameters, 
thereby varying Bond number, and photographing the meniscus shape. 

The measured meniscus was corrected for distortion by 
calculating correction factors based on basic laws of reflection and 
refraction. The calculated correction factors were checked and veri- 
fied by reading photographs of ball bearings with known dimensions. 

This procedure yielded accurate results except very near the cylinder 
walls where distortion was greatest. As mentioned in Reference 16, 
because of distortion problems, difficulties were encountered in deter- 
mining exactly where the interfaces intersected the container walls, and 
accurate contact angle measurements were not possible. 

Lockheed measured surface shapes for Bond numbers ranging 
from 8 to 53 using three test liquids: water, carbon tetrachloride, and 
methyl alcohol. Contact angles of 66„ 18, and 17 degrees were speci- 
fied for the water, carbon tetrachloride, and alcohol, respectively. 
However, the contact angles specified for carbon tetrachloride and 
alcohol are believed to be incorrect for the following reasons: 
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1. Lockheed investigators did not express confidence in the 
mea sirred contact angles. 

2. Contact angles of zero degrees are usually quoted in litera- 
ture for carbon tetrachloride and methyl alcohol in contact with glass 
or Lucite. (References 17 and 18.) 

3. The theoretical shapes calculated by Lockheed indicated 
that the actual contact angles were lower than the measured values. 

The measured contact angle for water should be correct since 
the measurement accuracy for such large angles should be good. 

Using contact angles of zero degrees for carbon tetrachloride 
and alcohol and 66 degrees for water, theoretical surface shapes were 
determined and are compared with the Lockheed experimental data in 
Figures 16 through 18. As illustrated in these figures, the theoretical 
profiles agree exceptionally well with the experimental data. In fact, 
if the actual contact angle is known, it is believed that the static equili- 
brium interface profiles can be computed with greater accuracy than 
they can be measured due to the distortion and reflection problems 
inherent in such experimental measurements. 

Attempts have been made to obtain additional experimental data 
in actual low gravity environments provided by the Marshall Space 
Flight Center (MSFC) drop tower facility. Cylindrical and spherical 
containers six inches in diameter were utilized with petroleum ether, 
a zero contact angle liquid, as the test fluid. However, when a liquid- 
vapor system is subjected to a sudden decrease in acceleration, such 
as that encountered in drop tower testing, certain interface oscillations 
must occur before the equilibrium configuration is attained. Theoreti- 
cal and experimental evaluations of such interface oscillations have 
been presented by Paynter, Fung, and Siegret, et al in References 19, 
20, and 21, respectively. 
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The maximum test time available in the MSFC drop tower 
(4 . 3 seconds) was insufficient to permit the interface to attain complete 
equilibrium. Thus, at best, the surface profiles attained only a state 
of quasi-equilibrium and the data can not be used to accurately verify 
theoretical interface solutions. However, based on preliminary 
comparisons, it can be stated that the experimental profiles did 
appear to oscillate about the theoretical static equilibrium shapes. 

It is anticipated that this drop tower data will be published in a MSFC 
document in the near future. 
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FIGURE 16. EXPERIMENTAL AND THEORETICAL 
INTERFACE SHAPES IN A CYLINDER 
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Dimensionless Height 


— Theoretical Zero Contact Angle Surface 
O Experimental Data From Ref. 16 
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FIGURE 17, EXPERIMENTAL AND THEORETICAL 
INTERFACE SHAPES IN A CYLINDER 
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Dimensionless Height 


— Theoretical Surface For Contact Angle = 66® 
O Experimental Data From Ref. 16 
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Test Liquid: Water 
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FIGURE 18. EXPERIMENTAL AND THEORETICAL 
INTERFACE SHAPES IN A CYLINDER 
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VII. CONCLUSIONS 


Based on the theoretical and experimental data acquired by these 
and other investigators, the following conclusions are made. 

1. Use of a polar coordinate system in developing the basic 
differential equation for liquid-vapor interface profiles eliminates con- 
vergence difficulties encountered in the solution of previously developed 
differential equations. 

2. The polar coordinate system enabled the incorporation of a 
Bond number based on a container dimension into the basic differential 
equation as opposed to a Bond number based on interface radius of 
curvature. 

3. Using the Runge-Kutta numerical technique, the interface 
equation developed herein can be readily solved by a computer. The 
method imposes no significant limitations on contact angle or Bond 
number . 

4. The effect of contact angle on surface shapes decreases with 
increasing Bond number and becomes negligible as the zero degree con- 
tact angle liquid surface approaches flatness. 

5. The theoretical equilibrium interface profiles, determined 
using the methods presented herein, correlate well with experimentally 
measured surface profiles. 

6. Due to the distortion and reflection problems inherent in 
measuring actual surface profiles and contact angles, it appears that 
the profiles can in most cases be theoretically calculated with greater 
accuracy than the surfaces can be measured. 
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PRECEDING PAGE BLANK NOT PlLMED* 


APPENDIX A 

LIQUID. VAPOR- SOLID INTERFACES 

That the shape of a droplet can be significantly affected by the 
presence of a solid is a well known observation. The degree to which 
the liquid- vapor interface shape is influenced is dependent on whether 
the cohesive* or adhesive* forces dominate, that is, the "wettability". 

The term "wettability", sometimes called "spreadability" is easily 
illustrated by noting a common, everyday occurrence; the effects of 
waxed and unwaxed surfaces on liquid droplet behavior. Droplets of 
water on the waxed surface will form "beads", while water on the 
unwaxed surface rapidly spreads or wets the surface. Whenever the 
degree of attraction between the liquid and solid is discussed (wetta- 
bility) the term usually invoked is "contact angle". 

Contact angle, as described in Figure 1A, is the angle (measured 
in the liquid) between the solid-liquid and the liquid- vapor interfaces. 

If a contact angle less than 90 degrees exists, the surface is said to 
be wetted; a contact angle greater than 90 degrees denotes a "non-wetting" 
of the surface. Many liquid- solid surfaces demonstrate total wetting, 
i.e., contact angles of zero degrees; but it is impossible to have a 
perfectly non- wetting liquid-solid surface, i.e., a contact angle equal 
to 180 degrees. The only liquid approaching complete non -wettability 
is mercury, which has a contact angle of about 125 degrees on glass 
surfaces. 


* These terms are defined in this study under "Basics of Liquid-Vapor- 
Solid Systems". 
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FIGURE 1A. CONTACT ANGLE MEASUREMENT 

Over 150 years ago, Thomas Young proposed treating the contact 
angle of a liquid as the result of the mechanical equilibrium of three 
surface tensions acting on the line of contact between a liquid- gas 
interface and a solid surface. These surface tensions are usually 
termed the solid-vapor (<r sv ), liquid- solid (a^ s ), and liquid-vapor (cr, g v ) 
surface tensions and are assumed to act in a direction parallel to each 
of their respective interfaces at the line of contact (see Figure 2A). 

This line of contact can be displaced to increase the solid-gas interface 
at the expense of the solid-liquid interface. If the solid-liquid surfaces 
exerted no force upon the line of contact, then obviously no equilibrium 
position would be possible since a force, cr^ v cos 0, parallel to the solid 
surface acts on this line. 

There must, therefore, be forces of the same nature as surface 
tensions that act through the line of contact, P, and are associated 
with the solid-vapor and solid-liquid interfaces. Young (Ref. 22) 
proposed that contact angle be related to the surface tensions of the 
three surfaces by the relation 


cr. cos a = 
i v 


sv 


- cr 


si 


(1A) 
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FIGURE 2A. SURFACE TENSIONS AT A LIQUID- VAPOR- 
SOLID INTERFACE 

Since the solid-vapor and solid- liquid surface tensions are not well 
understood and are difficult if not impossible quantities to measure 
experimentally. Young's equation is an extremely useful tool in that it 
expresses the effects of these surface tensions in terms of measurable 
quantities, that is, contact angle and liquid- vapor surface tension. 
However, the relation is deceptively simple and has been the 
source of many arguments. Biker man (Ref. 23) has criticized the 
equation on the grounds that the equilibrium conditions are discussed 
only with respect to forces parallel to the surface, and that no account 
is taken of the component 0-4 v sin a normal to the solid surface. 

A clear statement of the problem and a thermodynamic justifica- 
tion of Young's relation was given by Johnson (Ref. 24). Also, 

Lester (Ref. 25) has recently given a sophisticated treatment of the 
equation and showed that it is correct so long as the solid is not 
easily deformable. 

In addition, the surface condition of a solid can significantly 
affect the contact angle of a liquid on the solid. One such effect, 
which was analyzed by Wenzel (Ref. 26), is surface roughness. It 
seems that liquid on a rough surface will exhibit lower contact angle s than on 
a smooth one, because the surface irregularities provide many capillary 
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paths for the liquid and thereby cause the liquid to spread. Wenzel 
suggested a modified form of equation (1A) 

cr lv cos a = r(<r sv - <r s | ) (2A) 

where a is the average apparent contact angle and r is the ratio of true 
to apparent area of the solid. 

A second influence on contact angle is the presence of molecules 
adsorbed at the solid interface. This influence was observed by 
Langmuir (Ref. 27) when he measured a contact angle hysteresis 
effect whenever the line of contact between the liquid and solid was in 
motion, that is, the measured contact angle depended on whether the 
boundary was advancing or receding. Langmuir attributed this varia- 
tion in contact angle to a monomolecular layer that was adsorbed when 
liquid advanced over the surface and thereby decreased the contact 
angle when the liquid receded. 

In conclusion, the use of Young's equation and the concept of 
contact angle appears to be valid if the restrictions involved are 
properly understood. If further information on liquid- solid- vapor 
interfaces is desired. Reference 28 contains a comprehensive summary 
of recent works on the subject. 
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APPENDIX B 


BOUNDARY CONDITIONS DEPENDENT ON 
CONTAINER SHAPE 


The expressions for boundary conditions that are dependent on 
container shape are developed in this Appendix. As noted previously, 
the liquid- vapor interface shapes are symmetrical about the vertical 
axis, and therefore, the relations for container boundary and contact 

I 

angle can be formulated based on a vertical cross-section of the 
container . 

A. Prolate and Oblate Spheroids 

Referring to Figure IB, the equations for container boundary 
(Yg) contact angle (a), and empty fraction (|3) are derived as follows: 

1, Container Boundary - In this study, the equation for a vertical 
cross-section of both a prolate and oblate spheroid (an ellipse) is non- 
dimens ionalized with respect to the vertical semi-axis of the ellipse, 
i.e., the Z/R Q intercepts are always (0, -2). Therefore, the boundary 
for an ellipse is defined by the single relation 


Y - 2 b 2 cos 9 

B “ sin 2 0 + b 2 cos 2 0 

where "b" is always the horizontal dimension. Therefore, 


(IB) 


b < 1 Prolate Spheroid 
b > 1 Oblate Spheroid 

2. Contact Angle - As proven in Reference 29, the acute angle, 4*» 
between the tangent to the boundary, Yg = y(0), at point "P" and the 
line OP is defined by the relation 
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.g 


Bond No. 



FIGURE IB. GEOMETRY FOR PROLATE AND OBLATE 
SPHEROID BOUNDARY CONDITIONS 


tan 4* s - 

y'(e) 


(2B) 


Since 


/B , 2 b 2 cos Q 

^ ~ sin 2 0 + b* cos 2 0 


Performing the operation indicated in equation (2B) yields 

, f cos 9 (sin 2 9 + b 2 cos 2 8) 1 

y - arc n 8 lxi^ ©' — b 2 cos* 9 sin 0 + 2 sin0 cos 2 0j 


(3B) 


Also, note that since X * R 0 y sin 0 and Z a - R Q y cos 0 then the 
angle, y, measured from the vertical MP and the tangent to y(0) at "P" 
can be specified as follows 


y a arc tan 
y a arc tan 


(if) at p° int p 

/y 1 sin0 + y cos 0\ 
\y sin0 - y 1 cos 0/ 


(4B) 


From the geometry of Figure IB, it is apparent that the contact angle 
a is 


aay+<|>ay + 0- vJj 


(5B) 


Substituting in equation (5B) from relations (3B) and (4B) the final equa- 
tion for contact angle is 


a 


. arctan (? -■ ^ + .,, Y . <-° _ 3 8\ 
\y sin0 - y 1 cos 0/ 


+ 0 


* WWW V 1 W V | W WWW v I 

- arctan t | ' — , r 1 1 " • * ' . — : — • • . — rr 

_sxn 3 0 - b 2 cos 2 0 sm0 + 2 sin0 cos 2 0 


cos 0 (sin 2 0 + b 2 cos 2 0) 


(6B) 


3. Empty Fraction - Since the equation for an ellipse with the origin 
at the top of the ellipse is 

p- + (Z + l) 2 * 1 


59 



Then it is easily shown that for the total container volume (V t ) and 
vapor volume (V v ), the relations are 


V t = 2 it 


I 


X 1 2 dZ 


0 

J b 2 [1 - (Z + l) 2 ] dZ 
-1 


V t 



IT b 2 


and 


0 

Y v = tt J b 2 [ 1 - (Z + I) 2 ] dZ 



(7B) 


(8B) 


and empty fraction, (3, is simply 


P 


“v t 



(9B) 


B. Spheres 

It is a very simple matter to develop relations for the boundary 
conditions in a spherical container by setting the semimajor and semi- 
minor axes equal, i. e. , set b = 1, in the equations for the spheroidal 
containers to make the relations applicable to a spherical cross section 
(circle). On inserting b = 1 in relations (IB), (6B), and (9B), 


1. Container Boundary 
Yg = 2 cos 0 

2. Contact Angle 


^ „ . 180 , /y' sin 6 + y cos 0\ 

a = 2 0+ arctan l - 2 — ~) - 90 

it \y sin 0 - y* cos 0/ 


(10B) 


(HB) 
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3. Empty Fraction 

Since the equation for empty fraction in the spheroids is 
independent of "b", relation (9B) need not be changed for application 
to spherical containers. 

C. Cylinders 

The cylinder is the most simple geometrical configuration to 
analyze from the standpoint of interface shape since the container 
boundary is independent of the vertical coordinate "Z". Therefore, 
using Figure 2B to illustrate the geometry involved, it is easy to 
show that the boundary conditions for a cylinder are as follows: 

1. Container Boundary - The vertical cross section of a cylindrical 
container that is symmetric about the vertical axis and of infinite 
height is prescribed by 



(12B) 


2, Contact Angle - Contact angle is the slope of a tangent to the liquid- 
vapor inter face /container boundary or 

/dX\ 

a = arctan ("jT^n at point P 


and. 


a = arctan 


r y' sin 8 + y cos 9\ 
vy sin 9 - y 1 cos 9/ 


(13B) 


3. Original Liquid Height - Since the interface shape in a cylinder of 
at least one container radius in depth is independent of the empty frac- 
tion, a calculation of (3 is not necessary to specify the interface in a 
cylinder. However, in order to designate the high gravity liquid level 
corresponding to the low gravity interface, it is necessary to compute 
the liquid volume participating in the interface deformation. The equa- 
tion employed is 
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Bond No. = 00 



FIGURE 2B. GEOMETRY FOR BOUNDARY CONDITIONS 
IN A CYLINDER 



d< v l ) * 



0 X 

= J ir (1 - y 2 
0 


sin 2 0)(y sin0 - y' cos0)d0 


Therefore, the original liquid height "h" is 


h 



e l 

^ (1 - y 2 sin 2 8)(y sin 0 - y' cos 0) d0 

0 


(14B) 


4. Special Boundary Conditions Applicable to a Cylinder - Since the 
liquid fill level has no significance in the calculation of interface shape 
in cylinders, the only requirement for the observation point position 
(coordinate system origin) is that it be above the interface* Therefore, 
the distance to the low gravity face can be set equal to one, that is, 
y Q as 1, and the initial conditions for the main interface differential 
equation (equation (37)) can be simplified to 


y(o) = y Q = 1 

y'(0)=0 (15B) 

y” = 1 - K 0 
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APPENDIX C 


COMPUTER PROGRAM FOR DETERMINING 
INTERFACE SHAPE 


A. General Computer Procedure 

A general outline of the computer program sequence in engineering 
terms is as follows: 

1. Print out input data which includes: 


COMPUTER 

SYMBOL 

ENGINEERING 

SYMBOL 

DEFINITION 

UNITS 

KSVHC 


Container Designation: 

1 - Spherical 

2 - Vertical 

3 - Oblate Spheroid 

4 - Prolate Shperoid 


B 

b n 

Bond Number 

Dimensionless 

YG 

yg 

Distance from Origin to 
Surface for Bond Number = « 

Dimens ionle s s 

A2 


Spheroidal Container Major 
Axis 

Dimens ionles s 

B2 


Spheroidal Container Minor 
Axis 

Dimensionless 

CONANG 

a 

Contact Angle 

Degrees 

Y 0 * ** 

1 

' 

y o 

Estimated Distance from 
Origin to Low Gravity Sur- 
face Center Point 

Dimensionless 


2. The input data are used as initial values in the "Runge-Kutta" itera- 
tion solution of the main differential equation (33), which continues until the 

* The authors wish to acknowledge the contributions of Mrs Pam T. Hughes, 
of Computer Sciences Corporation, to the development of the computer 
program outlined herein. 

After some experience has been acquired in determining interface shapes, 
initial values of y Q can be determined that will enable more rapid computer 
convergence. Also, note that y Q = 1 in the case of a cylinder, since fill 
level has no influence on interface (see Appendix B). 
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interface intersects the container boundary. The Runge-Kutta formulas 
applied are those specified in Reference 11 for the integration of second- 
order equations of, the general type, 

y" =f(e, y, y') 
and consist of 


k l = A0 f(0 n> y n , y n ') 

k 2 = a© f(e n + y n + y n ’ + A 6 


L Ae Ae , a 0 , kA 

\ 0 n + T’ y n + T y n + T k * » y n + z) 

k 3 = A0 f(e n + y n + ~ y n ' + ^ k i> Yn + ^) 


A0 


k 4 = A0 f(e n + Ae, y n + Ae y n ' + — k 3 , y n ' + k 3 ) 


Ay = A0 + — (k 3 + k 2 + k 3 ) 


Ay' = (k x + 2k 2 + 2k 3 + k 4 ) 

3. Solve for interface /container boundary contact angle, a, and compare 
with the desired contact angle. If the calculated contact angle is negative, 
K Q is too large and must be decreased. However, if the angle is positive 
but too large, it is necessary to increase K Q . In either case, new values 
of K 0 are selected and used as new input in the Runge-Kutta solution until 
the desired contact angle is approached. 


4. Empty fraction* (vapor voluriie/ container volume), |3, is determined 
and compared with the desired empty fraction. If the correct (3 is not cal- 
culated, the value of y Q determined is entered in step 1, and new (3 is com- 
puted. This procedure is continued until enough data is generated to gener- 
ated to extrapolate or interpolate a curve fit of "yo versus (3" for the 
correct y Q . 

5. The correct value of y Q is entered in step 1 with the original value of 
K 0 , and the entire procedure is repeated until the desired contact angle 
and empty fraction are obtained. 


* In the case of a cylinder, the empty fraction criteria is by-passed. 
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6. Program output data is printed: 


COMPUTER 

SYMBOL 

ENGINEERING 

SYMBOL 

DEFINITION 

UNITS 

X/R; OX 

X/Ro 

Horizontal Distance fr om 
Container Vertical Axis 
to Liquid- Vapor Interface,- 

Dimen s ionl e s s 

Z/R; OZ 

Z/Ro 

Vertical Distance from 
X Axis to Liquid-Vapor 
Interface 

Dimensionless 

ANG; CONTACT 
ANGLE 


Contact Angle 

Degrees 

DXDY 

Ko 

Parameter Related to 
Curvature at Interface 
Centerpoint 

Dimensionless 

YO 

Yo 

Distance from Or igin 
to Low Gravity Surface 
Centerpoint 

Dimensionless 

TH; THETA 

0 

Angle Measured from 
Vertical Axis to y 

Degrees 

BETA; EMPTY 
FRACTION 


Calculated Empty 
Fraction 

Dimensionless 

BETAD 

D 

Desired Empty Fraction 

Dimensionless 

Y 

y 

Distance from Or igin to 
Liquid Surface 

Dimens ionles s 


B. Program Limitations 

As the low gravity interface shapes approach flatness at high Bond 
numbers, increasingly accurate values of K 0 are required because K 0 is 
approaching zero. In the present program, difficulty is encountered in 
obtaining contact angles of less than approximately five degrees at Bond 
numbers greater than or equal to 200, because the computer (GE 235) is 
unable to store the very small variations of K Q required for further con- 
vergence. This problem could be eliminated by using a more accurate 
computer (IBM 7094) if interface shapes at very high Bond numbers are 
desired; however, the interface shape variation caused by a five-degree 
contact angle deviation at Bond numbers above 150 is insignificant for 
most applications. 
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C. Programming Information 


1. Definition of Terms - Some of the significant terms not defined 
previously include: 



ENGINEERING 

SYMBOL 

DEFINITION ' 

UNITS 

DDY 

y" 

Second Derivative of y with 
Respect to 0 

Dimensionless 

DY 

y' 

First Derivative of y with 
Respect to 0 

Dimensionles s 

S 

Sin 0 

Sin 0 

Dimensionles s 

C 

Cos 0 

Cos 0 

Dimensionless 

VY 

V L 

Liquid Volume 

Dimens ionle s s 

DVB 

av l 

Incremental Liquid Volume 

Dimensionless 

DBET 

Ap 

Incremental Empty Fraction 

Dimensionless 
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2. Flow Chart 
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1 
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2 
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220 
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3. Program Listing 


C SPHERE INTERFACE* EMPIRICAL - VFRSION C-2 

KOMMON 0X12000) »OZ< 2000) 

COMMON ANSI 3) »ARG(3) »TRDX(2) » TBANG ( 2 ) 

29 FORMAT (12) 

30 FORMAT I5E16.8) 

31 FORMAT ( 6HX/R = »F7.4»8H» Z/R = *F7.4) 

10 FORMAT ( 5X » ) n HNO CONTACT) 

66 FORMAT ( 1 4HC0NT ACT ANGLF *F7.3»RH DPGRPES) 

34 FORMAT f //»1X3H3= El 1 .5 , 3X6HBETAD-EH . 5 » 3X3HYG=E1 1 . 5 » 3X3HA = FI 
13X3HB= El I • 5 »2X11HANG DESIRED » 2XF1 1 .5 ) 

53 FORMAT ( 1H1 ♦ 1X27HSPHERE INTERFACE* EMPIRICAL) 

63 FORMAT (17HEMPTY FRACTION = *F6.4) 

993 FORMAT ( 1H1 » 16HPROL ATE SPHEROID) 

994 FORMAT ( 1 HI* 15HOBLATF SPHEROID) 

995 FORMAT! 1H1 »8HCYL I ND3R) 

996'-, FORMAT ( /?MB=F1 1 ,5 ) 

D I=3«141 c; 92653S 
70 READ 37 » K SVHC 

READ 30,R,YG*A2,B2»C0NANG 
DTHD=360 • 

CDTH=.O03 

RETAD=( YG**?/4. ) *( 3.-YG) 

11=0 

IE(KSVHC-I )45,6,45 
6 <VHC=1 
PRINT 53 

PRINT 34 » B » BETAD » YG » A2 ,B2 *CON ANG 
GO TO 9 

45 KVHC=KSVHC-1 

GO TO(990, 991*992) »XVHC 

990 PRINT 99R 

PRINT 34»R»OFTAD»YG*A2*B? » CO NANG 
GO TO 9 

991 PRINT 994 

PRINT 34*B»BETAD»YG*A2»B2 , CONANG 
GO TO 9 

992 PRINT 995 
PRINT 996 »B 

9 DDDY=0. 

DFLX I N=1 • 

DDYO=.l 

DBC=.005 

YYO=YG 

37 FORMAT ( 1 1 ) 

I F ( SENSE SWITCH 19)501*500 
501 READ 30 » YYO » DDYO 
500 GO TO (215*216*271), KVHC 

215 Al=l * 

B 1 =B2 / A2 
GO TO 217 

216 Al=l. 

B1=A2/R2 

217 YO=YYO 
DYO=DDYO 

271 DELX=DELX I N 
DXDY=DDDY 

ddx=dxdy 

4999 JJ=0 


1*5* 
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49 DO 503 1=1,2000 
OX ( I ) = 0 • 

50^ oz( i )*n. 

GO TO (219,219,220) ,KVHC 
22° DDY=1.-DXDY 
Y0=1 • 

60 TO 221 

219 DDY=YO* ( l.-DXDY ) 

221 M=0 

GO TO (1000,1000,1001) ,KVHC 
1001 DTH=P I /DTHD 
GO TO 1002 
1000 DTH=P I /DTHD 
100? TH=0. 

T =0 

KCDTH=0 
DY-0 c 
BETA=0. 

FRFT = ° • 

RE=2 ( DDY-YO) /Y0**2-R*Y0 

Y=YO 

L=0 

vv-o . 

DVP-O, 

0 V A = 0 « 

672 1=1+1 

M-M+l 

C = COSF ( TH+DTH/2 • ) 

5 = SI MF ( TH+DTH/2 • ) 

CA=OTH*RRY 

Y A = Y +D TH*D Y / 2 . +DT H *C A / 8 . 

D Y A= D Y +C A / 2 * 

AA=YA**2 
AR=DYA**2 
AC=AA+AR 
AD=AC*SORTF (AC) 

CR = DTH* ( ( 2 • * A A+3 • *A P ) / YA-OYA*r* A(* / ( AA*F ) + ( R*YA*C+ n F ) *AR/YA ) 
YA = Y+DTH*DY/2'.+DTH*CA/8. 

D Y A = D Y +C q / 2 • 

AA=YA**2 
AR=DYA**2 
AC= AA+AB 
AD=AC*OQRTF( AC) 

CC-DTH*( ( 2.*AA+3.*AP ) / YA-DYA*C*AC/ ( A A*5 ) + ( R*Y A*C+R p ) *AD / Y A ) 
YA=Y+DTH*DY+DTH*Cr/2. 

DY A=DY +CC 
F=SINF( TH+DTH) 

C-COSF ( TH+DTH ) 

AA=YA**2 
AR=OYA**2 
AC= AA+AR 
AD=AC*OQRTF ( AC) 

CD=DTH* ( ( 2 •* A A+3 • *AP ) / YA-nYA*C* AC/ ( A A*9 ) + ( R*YA*C+ qp ) *AD/vA ) 

YZ=Y 

nYZ=OY 

THZ=TH 

OLDY =Y 

Y=Y+DTH* ( DY+ ( CA+CR+CC ) /6. ) 



DY=DY+ ( CA+? • * ( C3+CC ) +C0 ) /6 • 

IF(Y)673,673,682 
682 IF'(Y*.OLDY)65'>»65I»65.J 
65" IF(KCDTH1651 ,662,651 
65? DTH=CnTH 
KC0TH=1 
651 A A= Y*Y 

AR=0Y*OY 
AC=AA+AB 
AD=AC*SOPTF( AT) 

DDY= ( 2.*AA+?.*AB ) /Y-OY*C*AC/ ( AA*S ) + ( B*Y*C+BE ) *AO/Y 
IF (5FNSF SWITCH 18) 10? » in* 

103 PRINT 105»Y»DY»DDY,TH,CA,CB*CC»CD 
105 FORMAT (5X»8F12.5) 

104 TH=TH+OTH 
L = L+1 

GO TO (224*224*225)* KVHC 
225 IF(Y*S-1. )1 ,93,93 
224 CK1=7.*BT**2*C/(S**2+R1**2*C**? ) 

TF( Y-C<1 I 1,4,7 

1 IF(L-l) 18*18*? 

2 OX ( I )=Y*S 
OZ ( I )=Y*C 
L = 0 

18 GO TO (8,8,83 ) *KVHC 

83 DVR= ( 1 , — Y**2*S*S )*(Y*S-DY*C)*DTH /?. 

VV=VV+OVR +DVA 
OVA=nVR 
SZ = S 
GO TO 3 

8 DBET * ( • 75* ( Y*S-DY*C ) * ( Y*S ) **2 ) / ( B 1**2 ) 
BETA=BETA+0.5*(EBFT+DRFT)*DTH 
FBET=DBET 

sz=s 

R CZ = C 

IF ( TH* 180. /P 1-90.) 672, 673,67? 

673 PRINT 10 

DELX=DFLX/4. 

DXDY=DDX+DELX 
GO TO 49 

7 AA=YZ-2.*B1**2*CZ/ (SZ**2+B1**?*CZ**2) 
AB=Y-2.*B1**2*C/ ( S**2+P1 **2*C**2 ) 

THA= ( AB*THZ-AA*TH ) / ( AB—A A ) 

S=S INF ( THA ) 

C=COSF ( TH A ) 

Y=? .*B1**2*C/ ( S**2+BI **2*C**2 ) 
DY=DYZ+(DY-DYZ)*(THA-THZ)/(TH-THZ) 

TH=THA 

4 0X( 1+1 )=Y*S 
OZ ( 1+1 )=Y*C 
GO TO 91 

93 TE=TH-DTH*( Y*S-1. )/ ( Y*S-YZ*S? ) 

S=S I NF ( TF ) 

C=COSF(TF) 

OX ( 1+1 ) = Y*S 
0Z( 1+1 )=Y*C 
DTE=TF+DTH-TH 
VV=VV+DVA*DTF/DTH 



DYF=DYZ+DTE* ( DY-DYZ ) / DTH 
THA=TE 

91 ANG2=ANG1 
ANG1=AN6 

GO TO ( 1005 * 1005 * 1006 ) »KVHC 

1006 ANG=180./PI*(ATANF( (DYE*S+C/S)/(1.-DYE*C) ) ) 

GO TO 1007 

1005 YS I = ATANF( (2.*B1**2*C*(S**2+P1**2*C**2) )/(-(-7.*Rl**2*S**3+ 
] 2.*B1**4*C**2*S-4.*B1**2*S#C**2 ) 5 ) *180. /PI 

6AV!=]80.*TH/PI-YSI 

ALAM=180./PI*(ATANF( (Y*C+DY*S) / (Y*S-0Y*C> ) ) 

ang=gam+alam 

1007 IFi(ANG1+ANG2-2.*ANG)59»646»59 
646 IF(ANG)6'T»948.948 

59 IF (ANG) 60 ♦ 46 .5 1 

46 I F ( J J-l ) 6647 » 6646 > 66.47 

6647 I F ( SENSE SWITCH 15)6648.6649 

6648 PRINT 101 . ANG ♦ DXDY ♦ Y » THA 

6649 DDX=DXDY 
DELX=DELX/2. 

OXDY=DXDY— DELX 
GO TO 49 

6646 IF (ANG) 60 .849- .849 
849 I F ( ABSF ( ANG-CONANG ) -1 • ) 948 » 948 ♦ 949 
949 I F ( SENSE SWITCH 15)950.951 
95° PRINT 101 .ANG.nXDY.Y.THA 

951 TBDX ( 1 ) =DXDY 
TBDX ( 2 ) =DDX 
TBANG ( 1 ) = ANG 
TBANG ( 2 ) =ANG1 
NFRR=1 

CALL MLAG(2.2.C0NANG»TBANG.TBPX.TDX»NERR) 

DXDY=TDX 

I F ( SENSE SWITCH 15)952.49 

952 PRINT 2077. DXDY 

2077 FORMAT ( 1X4HDXDYE1 6 »8 ) 

GO TO 49 

60 IF (SENSE SWITCH 15) 61.62 

61 COST=2.*C 

PRINT 101 » ANG » DXDY . Y ♦ THA 

62 DFLX=DELX/4. 

DXDY = DDX + DFLX 
I KODE= 1 

GO TO 49 

51 IF ( ANG-CONANG) 46.9948 .9948 

9948 I F ( ABSF ( ANG-CONANG )-. 5 ) 948 . 948 ♦ 52 

52 COST=2»*C 

IF ( IKODE) 852.853.852 

852 DELXaDELX/4. 

I KODF = r ' 

853 CONTINUE 
JJ=1 

I F ( SENSE SWITCH 15)854.855 

854 PRINT 101 > ANG » DXDY »Y » THA 

101 FORMAT ( 1X3HANGE12 • 5 » 2X4HDXDYE1 6 • 8 » 2X1HYE1 2 • 5 , 2X2HTHE1 2 » 5 ) 

855 DDX = DXDY ’ 

DXDY = DXDY + DELX 
GO TO 49 
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948 GO TO { 850 *850*74 )»KVHC 
850 SS=S*S 

DBET=(.75*<Y*S-DY*C)*I Y*S)**2 )/ (Bl**2 J 
ZZZ=-OZ(I+l) 

BETA=BETA+0.5*(EBET+DBFT)*{TH-THZ)+(3«/4.*(ZZZ**3/3.+ZZZ**2) ) 

rr=BETAD-BETA 

Y01=Y0 

PRINT 667 »BETA » YO 

667 FORMAT ( 1X4HBETAE16 «8 » 3X2HYOE1 ? • 5 ) 

IF(ABSF(BB ) -DBC 174,74,47 
47 NERR=1 

I F (II) 32 * 33 » 32 
33 11=11+] 

ANSI 1 ) =YQ 
YO=YO+DYO 
ARG I 1 ) =BFTA 
DELX = DELX IN 
DXDY=DDDY 

ddx=dxdy 

GO TO 4999 
32 IF ( 11-1)35*36*35 
36 11=11+1 
ANS I 2 ) =YO 
ARG ( 2 ) =BFT A 
YO=YO+DYO 
DFLX=DFLX IN 
DXDY=DODY 

ddx=dxdy 

GO TO 4999 
35 I F ( I 1-2 ) 39 » 38 » 39 

38 ANSI 3 ) =YO 

I I = 1 T + l 

39 ARG I 3 ) =BFTA 

I F I ARG II) -ARG (3) ) 78 * 5 8 » 58 
58 AT=ARG 1 1 ) 

ARG I 1 ) =ARG I 3 ) 

ARG I 3 ) = AT 
TS= ANS 1 1 ) 

ANSI 1 ) = AMS I 3 ) 

ANSI 3 ) =TS 

78 I F I ARG 1 1 ) -ARG I 2 ) )87, 88*88 

88 PP=ARG I 1 ) 

ARG I 1 ) =ARG I 2 ) 

ARG I 2 ) =PP 
TS=ANS 1 1 ) 

ANSI 1 )=ANSI 2 ) 

ANSI 2)=TS 

87 I F I ARG I 2 ) -ARG I 3 ) ) 89 *90 *-9.0 
90 ST=ARG I 2 ) 

ARG I 2 ) -ARG I 3 ) 

ARG I 3 ) -ST 
T S=ANS I 2 ) 

ANSI 2 ) =ANS I 3 ) 

ANSI 3 ) =TS 

89 N E R R = 1 

CALL MLAG I 3 * 3 *BETAD* ARG* ANS » ANSE »NERR ) 

PRINT 9101, ANSE 

9101 FORMAT! IX 16HEXTRAP0LATED YO E12.5) 



IF ( ANSE ) 730 *730*731 

730 DYOsOYO+DDYO/2. 

DELX*DELX IN 
DXDY = DDDY 

do6<=dxdy 
11*0 
YO*YYO 
GO TO 4999 

731 YO=ANSE 

ARG{ 1 ) =ARG ( 2 ) 

ARG ( 2 ) =ARG ( 3 ) 

ANS ( 1 ) =ANS ( 2 ) 

ANSI 2 ) = AN$ ( 3 ) 

ANS(*)=YO 
DFLX = DFLXIN 
DXDYsDDDY 
DDX=DXDY 

100 FORMAT (1X2E16.8) 

GO TO 4999 
74 NN=M+1 

GO TO (732*732*733 )*KVHC 

732 00 81 I = 2 *M » 2 

81 PRINT 3 1 » OX ( I ) »0Z ( I ) 

PRINT 31* OX ( NN ) *OZ(NN) 

GO TO 734 

733 00 82 I* 2 * M *4 

82 PRINT 3 1 » OX ( I ) »0Z( I > 

PRINT 31 * OX ( NN ) ,OZ(NN) 

734 PRINT 66 » ANG 
THAA*THA#180./PI 
PRINT 117»0X0Y*Y0*THAA 
GO TO (870,870,860) *KVHC 

860 PRINT 1010»VV 
1010 FORMAT ( 2X7HHEIGHT»E1 2 • 5 ) 

GO TO 890 
870 PRINT 63 »BET A 

117 FORMAT ( 2X5HDXDY = F1 6 • 8 » 3HY0* El 6 • 8 »2X6HTHETA= E16.8) 
890 rALL ENOR 
GO TO 70 
END 
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D. Deck Setup 


1. Computer Configuration 


a. . 

Computer 

GE-235 

b. 

Core size required 

16 K 

c. 

Language 

F H 

d. 

System 

SLEM 

e. 

Plotter required 

NONE 

f . 

Card punch required 

NONE 


2. Estimated Run Time - Five minutes 

3. Restart Procedure - Multiple Cases Read by ENDR 

4. Deck Sequence 

ZERO MEMORY 16 K 
SEEM CALL CARD 
ID CARD 

BINARY FOR MAIN DECK 
BINARY FOR MLAG 
BINARY FOR ENDR 
1 /2 PUNCH 
DATA INPUT DECK 

5. Diagonistics 

If the program seems to be cycling, the program can be run with sense 
switch 15 down and the computer will print out the calculated contact 
angle, DXDY, and YO after each pass through the Runge-Kutta 
procedure. 
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6 . Input Sheet 


Job No. 

Card 1 
Card 2 


INPUT DATA SHEET 
Date : ... 

574260 

KSYHC = 1, 2, 3, 4 

1 - Spherical Shape, 2 - Vertical ellipse, 
3 - Horizontal ellipse, 4 - Cylinder 


Bond No* 

YG - Fill Level 

A 2 - Major Axis 

B 2 - Minor Axis 

Contact Angle 






Col 1 16 17 32 33 48 49 64 65 80 


Y 0 

DY 0 




Card 3 
OPTIONAL 


Col 1 


16 


17 


32 




7. Example Output - Sphere 


SPHERE INTERFACE# EMPIRICAL 


g« 0,140006 03 BlTAOiQ »58883E 00 YQaO,11198i 01 A« 0,100006 01 

NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 

NO contact 

NO CONTACT 
NO CONTACT 

SETA 0,578803226 00: yQ 0,111906 01 

NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 

BETA 0,032698466 00 TO 0*121906 01 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 

BETA 0,723317326 00 YO 0# 131988 01 


EXTRAPOLATED TO 
NO CONTACT 
NO CONTACT 
NO CONTACT 
N,Q CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 
NO CONTACT 


0 , 11323 * 01 


BETA 

0,988809706 

00 

YO 0,1 

X/R « 

0.0198, 

Z/R 


1.1323 

X/R ■ 

0 , 0399, 

Z/R 


1.1323 

X/R ■ 

0.0993, 

Z/R 


1,1323 

X/R ■ 

0,0792, 

Z/R 


1.1323 

X/R ■ 

0.0991, 

Z/R 


1,1323 

X/R « 

0.1190, 

Z/R 


1,1323 

X/R ■ 

0.1390, 

Z/R 


1,1323 

X/R ■ 

0.1591, 

Z/R 


1,1323 

X/R ■ 

0*1793, 

Z/R 


1,1323 

X/R ■ 

0.1997, 

Z/R 


1,1323 

X/R « 

0.2201# 

Z/R 


1.132J 

X/R ■ 

0.2407, 

Z/R 


1,1323 

X/R ■ 

0.2614, 

Z/R 


1,1323 

X/R * 

0.2823, 

Z/R 


1,1323 

X/R « 

0.3034, 

Z/R 


1,1323 

X/R « 

0.3247, 

Z/R 


1,1322 

X/R « 

0.3462, 

Z/R 


1,1322 

X/R « 

0.3679, 

Z/R 


1,1322 

X/R » 

0.3898, 

Z/R 


1,1322 

X/R « 

0.4121, 

Z/R 


1,1322 

X/R ■ 

0.4346, 

Z/R 


1,1321 

X/R * 

0,4374, 

Z/R 


1,1321 

X/R « 

0,4803, 

Z/R 


1,1321 

X/R * 

0.3040, 

Z/R 


1.1320 

X/R ■ 

0.5278, 

Z/R 


1,1319 

X/R « 

0.5520# 

Z/R 


1,1318 

X/R « 

0,3766, 

Z/R 


1,1316 

X/R ■ 

0.6016, 

Z/R 


1,1314 

X/R « 

0 . 6270# 

Z/R 


i,l3ll 

X/R « 

0.6526, 

Z/R 


1.1307 

X/R * 

0*6791, 

Z/R 


1,1302 

X/R * 

0.7058# 

Z/R 


1,1295 

X/R » 

0.7328# 

I/R 


1.1285 

X/R * 

0.76Q2, 

Z/R 


1*1271 

X/R * 

0.7879, 

Z/R 


1*1252 

X/R • 

0.8156# 

Z/R 


1,1226 

X/R * 

0.8433# 

Z/R 


1*1191 

X/R » 

0*3706# 

l/R 


1*1143 

X/R » 

0,8971# 

Z/R 


1*1079 

X/R « 

0.9225# 

Z/R 


1,0994 

X/R « 

0.9462# 

l/R 


1,0884 

X/R * 

0.9673# 

Z/R 


1,0745 

X/R * 

0.9856# 

Z/R 


1,0869 

X/R « 

0.9849# 

Z/R 


1,0556 

X/R • 

0.9896# 

l/R 


1,0480 

X/R « 

0*9938# 

Z/R 


1,0398 

X/R « 

0.9971# 

l/R 


1,0308 

X/R ■ 

0.9993# 

Z/R 


1,0208 

X/R ■ 

0 . • 

Z/R 


0 , 

X/R « 

0.9999# 

Z/R 


1,0128 

CONTACT ANCLE 

0*' 

0 

DECREES 

8X0Y» 0*3762990: 

■ 

03yo» 0, 


EMPTY FRACTION * 0*3888 


01 THETA* 0,446319966 02 


l« 0 .100008 01 ANQ 06SSRBD 0, 






PW 


PAGE BLANK 


hqt 
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